This notebook can be called with knit_notebook8_byESM.R; it can be called in a loop over ESMs and target experiments, and will knit an analysis html for each.
These values must be defined globally in knit_notebook8_byESM.R or uncommented and defined here to run this notebook directly:
# # what we will emulate:
# esm_name <- "MPI-ESM1-2-LR"
# esm_experiment <- 'ssp245'
#
# # How many times do we want to draw full generated ensembles:
# Ndraws <- 1
print(esm_name)
## [1] "MIROC6"
print(esm_experiment)
## [1] "ssp370"
We will only explore tol=0.1 degC for matching in this notebook (per CT, from the SD of the smoothed Tgav time series. Should check for each ESM, for now just use 0.1 value from CanESM).
# Then load the functions we will want to use, these are currently written in R and will be translated into python.
source("functions/nearest_neighbor_matching.R") # the function match_neighborhood is defined here
source("functions/permuting_matches.R") # the function permute_stitching_recipe is defined here
source("functions/stitching_functions.R") # the function stitch_global_mean is definend here
source("functions/gridded_recipe.R") # wrappers for the work flow from matching to output csv needed by python for
# layering in gridded data to output as netcdf
set.seed(1)
Load the inputs that will be used. In addition to SSP534-over having enough data issues to fully exclude, at least some of the realizations for SSP434 and SSP460 just don’t have future data at all for CanESM5. And for other models, some realizations have odd windows (Like a 2012-2014 window) and those realizations are excluded as well. We need to go back to pangeo and fix this for sure, but for now, I’ll work with the archive I have. Because we don’t want to lose too much of the archive now, I will only remove those realizations and not a full experiment for a single bad data point.
# Time series of the raw global mean temperature anomaly, this is the data that is
# going to be stitched together. Historical data is pasted into every scenario. No
# smoothing has been done.
tgav_data <- read.csv("inputs/main_raw_pasted_tgav_anomaly_all_pangeo_list_models.csv",
stringsAsFactors = FALSE) %>% dplyr::select(-X)
# A chunked smoothed tgav anomaly archive of data. This was tgav_data but several
# steps were taken in python to transform it into the "chunked" data and save it
# so that we do not have to repeate this process so many times.
read.csv(paste0('inputs/', esm_name, '_archive_data.csv'), stringsAsFactors = FALSE) %>%
# Exclude ssp534-over:
dplyr::filter(experiment != 'ssp534-over') %>%
# eliminate any individual realizationXexperiment runs that don't have all the years:
group_by(experiment, variable, model , ensemble) %>%
mutate(nWindows = n()) %>%
ungroup %>%
filter(nWindows == 28) %>%
select(-nWindows) ->
archive_data
# The target data - all ssp245 realizations (smoothed because it's the target for
# matching)
archive_data %>%
filter(experiment == esm_experiment) ->
target_data
# The corresponding unsmoothed original data for validation
tgav_data %>%
filter(model == unique(target_data$model) & experiment == unique(target_data$experiment)) %>%
filter(ensemble %in% unique(target_data$ensemble)) %>%
select(-file, -timestep, -grid_type) ->
original_data
We also need to prepare the file of path metadata for pangeo. This is a csv file that contains the pangeo paths to different variables of interest. By selecting the paths for tas, psl, and pr, keeping only those experimentXensemble combinations that have defined paths to netcdfs for all three variables, and then using the resulting experiment X ensemble combinations as an additional filter on the archives we use for matching, we make sure that we can then stitch together gridded data for all 3 of these variables no matter what our generated ensembles are. Adding or removing variables of interest is relatively trivial. For now, the path metadata file has information for tas, psl, pr, tasmin, and tasmax.
# The main pangeo file of path metadata so that we can filter the archives to include only
# the experiment*ensemble members with netcdfs for all variables of interest
read.csv("inputs/pangeo_path_metadata_tas_psl_pr_tmax_tmin.csv",
stringsAsFactors = FALSE,
na.strings = c("", " ", "NA")) %>%
rename(model = source_id, experiment = experiment_id, ensemble = member_id) %>%
filter(model == unique(target_data$model) ) %>%
select(model, experiment, ensemble, table_id, tas_zstore, psl_zstore, pr_zstore) %>%
# cut out any rows that don't have paths to netcdfs for all variables:
na.omit %>%
select(model, experiment, ensemble) %>%
unite("acceptable", c(model, experiment, ensemble), sep="~")->
acceptable_archive_entries_from_path_metadata
We’re actually going to generate multiple different ensembles for a target of ssp245, based on two different archives.
Archive 1 will exclude the ssp245 runs from potential matches.
Archive 2 will include the ssp245 runs - when we get to the point of targeting SSP585 runs, I think we will have to to include those runs in the archive to get at all decent generated Tgavs at the end of the century when temperature is higher. It’s still not just regurgitating the target ensembles, it will include mixing and matching from those parts.
We will also do some light exploration of the matching neighborhood size (tol argument).
# Archive data not including the target experiment:
archive_data %>%
# Only match to the same ESM:
dplyr::filter(model == unique(target_data$model)) %>%
# Don't match to any ensemble members in the target experiment:
dplyr::filter(experiment != unique(target_data$experiment) ) %>%
# Final filter to make sure we only use the experiment*ensemble members with pangeo
# path metadata for netcdfs for all varaibles of interest
unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
select(-acceptable) ->
archive_subset1
# what does this archive look like:
archive_subset1 %>%
select(experiment, ensemble) %>%
distinct() %>%
group_by(experiment) %>%
summarize(n_complete_ensemble_members = n()) %>%
ungroup %>%
kable()
| experiment | n_complete_ensemble_members |
|---|---|
| ssp119 | 1 |
| ssp126 | 50 |
| ssp245 | 3 |
| ssp434 | 1 |
| ssp460 | 1 |
| ssp585 | 50 |
# Archive data including the target experiment:
archive_data %>%
# Only match to the same ESM:
dplyr::filter(model == unique(target_data$model)) %>%
# Final filter to make sure we only use the experiment*ensemble members with pangeo
# path metadata for netcdfs for all varaibles of interest
unite("acceptable", c(model, experiment, ensemble), sep="~", remove=FALSE) %>%
filter(acceptable %in% unique(acceptable_archive_entries_from_path_metadata$acceptable)) %>%
select(-acceptable) ->
archive_subset2
# what does this archive look like:
archive_subset2 %>%
select(experiment, ensemble) %>%
distinct() %>%
group_by(experiment) %>%
summarize(n_complete_ensemble_members = n()) %>%
ungroup %>%
kable()
| experiment | n_complete_ensemble_members |
|---|---|
| ssp119 | 1 |
| ssp126 | 50 |
| ssp245 | 3 |
| ssp370 | 3 |
| ssp434 | 1 |
| ssp460 | 1 |
| ssp585 | 50 |
for (draw in 1:Ndraws){
# The matching to the archive without ssp245 results
match_neighborhood(target_data = target_data,
archive_data = archive_subset1,
tol=0.1) -> matched_data
matched_data %>%
# Convert them to a sample of individual Tgav Recipes:
permute_stitching_recipes(N_matches = 2000, matched_data = .,
archive = archive_subset1) ->
recipes_1_0p1
# convert these recipes to something that can be ingested by python
# to layer in gridded netcdfs and save
gridded_recipes_1_0p1 <- generate_gridded_recipe(recipes_1_0p1)
write.csv(gridded_recipes_1_0p1,
paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_not_in_matching.csv'),
row.names = FALSE)
# Proceed to stitch the actual Tgavs for analysis here.
recipes_1_0p1 %>%
# And stitch the recipes into tgavs:
stitch_global_mean(recipes=., data = tgav_data) %>%
# add some identifying info
mutate(tol = 0.1,
archive = paste0('without_', esm_experiment)) ->
rslts_1_0p1
# Check for collpase in these results
rslts_1_0p1 %>%
select(year, value, stitching_id) %>%
spread(year, value) ->
wide
# Collapse would have occured if any of these trajectories go through the same value
# in the same year -> if the length of distinct values in any column is less than
# the number of stitching_id's
counts <- apply(wide, 2, length)
print(paste('Archive without', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
# The matching to the archive with target experiment included results
match_neighborhood(target_data = target_data,
archive_data = archive_subset2,
tol=0.1) %>%
# Convert them to a sample of individual Tgav Recipes:
permute_stitching_recipes(N_matches = 2000, matched_data = .,
archive = archive_subset2) ->
recipes_2_0p1
# convert these recipes to something that can be ingested by python
# to layer in gridded netcdfs and save
gridded_recipes_2_0p1 <- generate_gridded_recipe(recipes_2_0p1)
write.csv(gridded_recipes_2_0p1,
paste0('generated_ensembles/', esm_name, '/gridded_recipes_for_python_', esm_experiment,'_draw', draw, '_experiment_in_matching.csv'),
row.names = FALSE)
recipes_2_0p1 %>%
# And stitch the recipes into tgavs:
stitch_global_mean(recipes=., data = tgav_data) %>%
# add some identifying info
mutate(tol = 0.1,
archive = paste0('with_', esm_experiment)) ->
rslts_2_0p1
# Check for collpase in these results
rslts_2_0p1 %>%
select(year, value, stitching_id) %>%
spread(year, value) ->
wide
# Collapse would have occured if any of these trajectories go through the same value
# in the same year -> if the length of distinct values in any column is less than
# the number of stitching_id's
counts <- apply(wide, 2, length)
print(paste('Archive with', esm_experiment,'; tol = 0.1: Is there any collapse in the generated ensemble?', any(counts != length(unique(wide$stitching_id)))))
# save it all off
write.csv(bind_rows(rslts_1_0p1,
rslts_2_0p1) %>%
mutate(draw = draw),
paste0('generated_ensembles/', esm_name, '/generated_ensembles_', esm_experiment,'_draw', draw, '.csv'),
row.names = FALSE)
}
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive without ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
## You have requested more recipes than possible for at least one target trajectories, returning what can
## [1] "Archive with ssp370 ; tol = 0.1: Is there any collapse in the generated ensemble? FALSE"
Comparing the generated ensembles to the training ensemble with error metrics along the lines of Tebaldi et al 2020.
That is, we need the means and variances of both the target ensemble and the generated ensemble. Specifically, mean across time and ensemble members, variance across time and ensemble members.
# values of the target data
mean_target <- mean(original_data$value)
mean_target_hist <- mean(filter(original_data, year <= 2014)$value)
mean_target_future <- mean(filter(original_data, year > 2014)$value)
var_target <- sd(original_data$value)
var_target_hist <- sd(filter(original_data, year <= 2014)$value)
var_target_future <- sd(filter(original_data, year > 2014)$value)
# generated ensemble files
# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name),
pattern = paste0('generated_ensembles_', esm_experiment),
full.names = TRUE)
for (i in 1:length(files)){
# Metrics for this draw #################################################################
generated_ensemble_draw <- read.csv(files[i], stringsAsFactors = FALSE)
# Values of the generated data and then the metrics
generated_ensemble_draw %>%
group_by(tol, archive, draw) %>%
summarize(mean_gen = mean(value),
var_gen = sd(value)) %>%
ungroup %>%
mutate(time_period = 'hist+future',
mean_target = mean_target,
var_target = var_target,
bias = abs(mean_target - mean_gen),
E1 = bias/var_target,
E2 = var_gen/var_target) ->
metrics
generated_ensemble_draw %>%
filter(year <= 2014) %>%
group_by(tol, archive, draw) %>%
summarize(mean_gen = mean(value),
var_gen = sd(value)) %>%
ungroup %>%
mutate(time_period = 'hist',
mean_target = mean_target_hist,
var_target = var_target_hist,
bias = abs(mean_target - mean_gen),
E1 = bias/var_target,
E2 = var_gen/var_target) ->
metrics_hist
generated_ensemble_draw%>%
filter(year > 2014) %>%
group_by(tol, archive, draw) %>%
summarize(mean_gen = mean(value),
var_gen = sd(value)) %>%
ungroup %>%
mutate(time_period = 'future',
mean_target = mean_target_future,
var_target = var_target_future,
bias = abs(mean_target - mean_gen),
E1 = bias/var_target,
E2 = var_gen/var_target) ->
metrics_future
generated_ensemble_draw %>%
select(stitching_id, tol, archive, draw) %>%
distinct %>%
group_by(tol, archive, draw) %>%
summarise(n_stitched = n()) %>%
ungroup ->
tmp
bind_rows(tmp %>%
left_join(metrics, by = c('tol', 'archive', 'draw')),
tmp %>%
left_join(metrics_hist, by = c('tol', 'archive', 'draw')) ,
tmp %>%
left_join(metrics_future, by = c('tol', 'archive', 'draw'))) ->
plotting_tbl
rm(tmp)
write.csv(plotting_tbl %>%
mutate(draw = draw),
paste0('generated_ensembles/', esm_name ,'/metrics_', esm_experiment,'_draw', i, '.csv'),
row.names = FALSE)
}
| tol | archive | draw | n_stitched | mean_gen | var_gen | time_period | mean_target | var_target | bias | E1 | E2 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.1 | with_ssp370 | 1 | 10 | 0.1441791 | 0.9791005 | hist+future | 0.1463943 | 0.9853313 | 0.0022152 | 0.0022482 | 0.9936764 |
| 0.1 | without_ssp370 | 1 | 30 | 0.1502757 | 0.9795846 | hist+future | 0.1463943 | 0.9853313 | 0.0038814 | 0.0039392 | 0.9941677 |
| 0.1 | with_ssp370 | 1 | 10 | -0.4648884 | 0.2535008 | hist | -0.4652324 | 0.2543384 | 0.0003440 | 0.0013527 | 0.9967064 |
| 0.1 | without_ssp370 | 1 | 30 | -0.4603009 | 0.2504051 | hist | -0.4652324 | 0.2543384 | 0.0049316 | 0.0193897 | 0.9845350 |
| 0.1 | with_ssp370 | 1 | 10 | 1.3127389 | 0.7726450 | future | 1.3198643 | 0.7833302 | 0.0071254 | 0.0090963 | 0.9863592 |
| 0.1 | without_ssp370 | 1 | 30 | 1.3217308 | 0.7698456 | future | 1.3198643 | 0.7833302 | 0.0018665 | 0.0023827 | 0.9827856 |
In general, all of our emulations have E1 close to 0 and E2 close to 1. As always, the question is sort of ‘how close is close enough’.
Again, it’s good to remember tgat because the draw of generated ensemble Tgavs is random, there’s no guarantee that for an individual draw, using the archive including the SSP245 points will guarantee better fits. It will certainly create a larger potential space of better fits to draw from. For SSP585, I suspect this will be more of a factor.
| tol | archive | draw | n_stitched | mean_gen | var_gen | time_period | mean_target | var_target | bias | E1 | E2 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.1 | with_ssp370 | 10 | 9 | 0.1480769 | 0.9752615 | hist+future | 0.1463943 | 0.9853313 | 0.0016825 | 0.0017076 | 0.9897803 |
| 0.1 | without_ssp370 | 10 | 30 | 0.1488753 | 0.9769803 | hist+future | 0.1463943 | 0.9853313 | 0.0024810 | 0.0025179 | 0.9915246 |
| 0.1 | with_ssp370 | 10 | 9 | -0.4590737 | 0.2483653 | hist | -0.4652324 | 0.2543384 | 0.0061587 | 0.0242146 | 0.9765149 |
| 0.1 | without_ssp370 | 10 | 30 | -0.4587786 | 0.2506305 | hist | -0.4652324 | 0.2543384 | 0.0064538 | 0.0253750 | 0.9854213 |
| 0.1 | with_ssp370 | 10 | 9 | 1.3129589 | 0.7700946 | future | 1.3198643 | 0.7833302 | 0.0069055 | 0.0088155 | 0.9831035 |
| 0.1 | without_ssp370 | 10 | 30 | 1.3147230 | 0.7729910 | future | 1.3198643 | 0.7833302 | 0.0051414 | 0.0065635 | 0.9868010 |
We want to know the average errors across many draws so that we can get a sense of the error of our method and how that relates to hyperperameters:
# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
pattern = paste0('metrics_', esm_experiment),
full.names = TRUE)
generated_metrics <- list()
for (i in 1:length(files)){
generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)
generated_metrics %>%
group_by(tol, archive, time_period) %>%
summarize(aggregate_E1 = mean(E1),
aggregate_E2 = mean(E2)) %>%
ungroup %>%
mutate(total_draws = length(unique(generated_metrics$draw))) ->
aggregate_metrics
kable(aggregate_metrics)
| tol | archive | time_period | aggregate_E1 | aggregate_E2 | total_draws |
|---|---|---|---|---|---|
| 0.1 | with_ssp370 | future | 0.0100603 | 0.9869164 | 10 |
| 0.1 | with_ssp370 | hist | 0.0131077 | 0.9920320 | 10 |
| 0.1 | with_ssp370 | hist+future | 0.0030050 | 0.9928589 | 10 |
| 0.1 | without_ssp370 | future | 0.0024767 | 0.9828202 | 10 |
| 0.1 | without_ssp370 | hist | 0.0263855 | 0.9856265 | 10 |
| 0.1 | without_ssp370 | hist+future | 0.0040445 | 0.9921914 | 10 |
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
geom_point() +
facet_wrap(~time_period, nrow =3) +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 0)
Pull off the years that end each window and that start each window; these are going to be years we pay particular attention to.
window_start_years <- unique(archive_subset1$start_yr)
window_end_years <- unique(archive_subset1$end_yr)
We are particularly concerned about the ‘jumps’ at each of these years being ‘wrong’. And the idea is that they would be wrong relative to both the target data and, I think, the generated ensembles in between the jump years.
So let’s work with the jumps themselves for validating. Take the year to year diff within each realization (target or generated ensemble member). We’ll look at the diff and abs(diff). We’ll visualize the jumps with a scatter view and then get into more quantitative metrics. I spent a while going back and forth between starting from fundamentals like this and tracking down a package that looks at more sophisticated methods. As I started reading through packages, I ended up deciding to stick to fundamentals because they’ll translate from R to python better and will provide a more solid understanding of what’s going on before we look at anything more sophisticated.
# read in all generated ensembles.
files <- list.files(path = paste0('generated_ensembles/', esm_name),
pattern = paste0('generated_ensembles_', esm_experiment),
full.names = TRUE)
generated_tgavs <- list()
for (i in 1:length(files)){
generated_tgavs[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_tgavs <- do.call(rbind, generated_tgavs)
# make a helper function for calculating all the jumps because you have Nyears-1 jumps and
# dply functions don't like that.
get_jumps <- function(data){
data.frame(left_year = data$year[1:(length(data$year)-1)],
jump = diff(data$value),
abs_jump = abs(diff(data$value))) %>%
mutate(variable = unique(data$variable),
stitching_id = unique(data$stitching_id),
tol = unique(data$tol),
archive = unique(data$archive),
draw = unique(data$draw),
jump_year = if_else(left_year %in% window_end_years, 'window_transition', 'in_window')) ->
x
return(x)
}
# calculate the yearly jumps
generated_list <- split(generated_tgavs,
list(generated_tgavs$tol,
generated_tgavs$archive,
generated_tgavs$draw,
generated_tgavs$stitching_id),
drop=TRUE)
lapply(generated_list, get_jumps) %>%
do.call(rbind, .) ->
generated_jumps
target_list <- split(original_data,
list(original_data$model,
original_data$experiment,
original_data$ensemble),
drop = TRUE)
lapply(target_list, get_jumps) %>%
do.call(rbind, .) ->
target_jumps
Generated jump year values during window transitions vs all target jump year values:
ggplot() +
geom_point(target_jumps, mapping = aes(x = left_year, y = jump), alpha = 0.5) +
geom_point(filter(generated_jumps,
jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.1, color='red') +
facet_grid(tol~archive)
ggplot() +
geom_point(filter(target_jumps,
jump_year == 'window_transition'), mapping = aes(x = left_year, y = jump), alpha = 0.5) +
geom_point(filter(generated_jumps,
jump_year == 'window_transition'), mapping = aes(x = left_year-1, y = jump), alpha = 0.1, color='red' ) +
facet_grid(archive~tol) +
ggtitle('note the years are staggered just for visual clarity')
What’s the appropriate comparison though?
From the perspective of the target data, there’s nothing special about those transition years, which basically knocks out 3. I’ll plot it for completeness but I definitely don’t think it’s the right comparison
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years, full ensemble (target or generated), \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nTotal number of draws for generated ensembles ', length(unique(generated_jumps$draw))))
focus_draw <- 1
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))
focus_draw <- 2
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years, full ensemble (target or generated), \nFocusing on a single draw of generated ensemble ', focus_draw))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('All years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))
ggplot() +
geom_density(target_jumps %>% mutate(id= 'target')%>% filter(jump_year == 'window_transition'), mapping = aes(x = jump, fill = id), alpha =0.2) +
geom_density(generated_jumps %>% mutate(id = 'generated') %>% filter(jump_year == 'window_transition', draw == focus_draw), mapping = aes(x = jump, fill = id), alpha = 0.2) +
facet_grid(archive~tol) +
ggtitle(paste('Window years target ensemble, \nWindow years generated ensemble, \nFocusing on a single draw of generated ensemble ', focus_draw))
Effectively quantifying how similar vs not the above distributions are, since they’re all pretty normal looking. I suppose we could do a KS test to support further if needed?
So we will compare all the years of generated data to all the years of target data, and the window transition years of generated data to all the years of target data. Individual draw and aggregate across draws errors
Individual draws:
# values of the target data
mean_target <- mean(target_jumps$jump)
var_target <- sd(target_jumps$jump)
for (drawID in 1:Ndraws){
# Metrics for this draw #################################################################
generated_jumps %>%
filter(draw == drawID) ->
generated_ensemble_draw
# Values of the generated data and then the metrics
generated_ensemble_draw %>%
group_by(tol, archive, draw) %>%
summarize(mean_gen = mean(jump),
var_gen = sd(jump)) %>%
ungroup %>%
mutate(mean_target = mean_target,
var_target = var_target,
bias = abs(mean_target - mean_gen),
E1 = bias/var_target,
E2 = var_gen/var_target) ->
metrics
generated_ensemble_draw %>%
select(stitching_id, tol, archive, draw) %>%
distinct %>%
group_by(tol, archive, draw) %>%
summarise(n_stitched = n()) %>%
ungroup %>%
left_join(metrics, by = c('tol', 'archive', 'draw')) ->
plotting_tbl
write.csv(plotting_tbl %>%
mutate(draw = drawID),
paste0('generated_ensembles/', esm_name ,'/jump_errors_allyears_', esm_experiment, '_draw', drawID, '.csv'),
row.names = FALSE)
}
# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/',esm_name , '/jump_errors_allyears_',esm_experiment,'_draw1.csv'),
stringsAsFactors = FALSE)
kable(plotting_tbl)
| tol | archive | draw | n_stitched | mean_gen | var_gen | mean_target | var_target | bias | E1 | E2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.1 | with_ssp370 | 1 | 10 | 0.0137593 | 0.1830368 | 0.0147099 | 0.1754115 | 0.0009507 | 0.0054196 | 1.043471 |
| 0.1 | without_ssp370 | 1 | 30 | 0.0134830 | 0.1819347 | 0.0147099 | 0.1754115 | 0.0012270 | 0.0069948 | 1.037188 |
ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
geom_point() +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 0) +
ggtitle('Error metrics for single Draw 1, comparing all \nyears of generated jumps to all years of target jumps')
Aggregated across draws:
# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
pattern = paste0('jump_errors_allyears_', esm_experiment),
full.names = TRUE)
generated_metrics <- list()
for (i in 1:length(files)){
generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)
generated_metrics %>%
group_by(tol, archive) %>%
summarize(aggregate_E1 = mean(E1),
aggregate_E2 = mean(E2)) %>%
ungroup %>%
mutate(total_draws = length(unique(generated_metrics$draw))) ->
aggregate_metrics
kable(aggregate_metrics)
| tol | archive | aggregate_E1 | aggregate_E2 | total_draws |
|---|---|---|---|---|
| 0.1 | with_ssp370 | 0.0048429 | 1.048425 | 10 |
| 0.1 | without_ssp370 | 0.0075740 | 1.040463 | 10 |
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
geom_point() +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 0) +
ggtitle(paste('Aggregate errors (jumps for all years), aggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))
Individual draw:
# values of the target data
mean_target <- mean(target_jumps$jump)
var_target <- sd(target_jumps$jump)
for (drawID in 1:Ndraws){
# Metrics for this draw #################################################################
generated_jumps %>%
filter(draw == drawID,
jump_year=='window_transition') ->
generated_ensemble_draw
# Values of the generated data and then the metrics
generated_ensemble_draw %>%
group_by(tol, archive, draw) %>%
summarize(mean_gen = mean(jump),
var_gen = sd(jump)) %>%
ungroup %>%
mutate(mean_target = mean_target,
var_target = var_target,
bias = abs(mean_target - mean_gen),
E1 = bias/var_target,
E2 = var_gen/var_target) ->
metrics
generated_ensemble_draw %>%
select(stitching_id, tol, archive, draw) %>%
distinct %>%
group_by(tol, archive, draw) %>%
summarise(n_stitched = n()) %>%
ungroup %>%
left_join(metrics, by = c('tol', 'archive', 'draw')) ->
plotting_tbl
write.csv(plotting_tbl %>%
mutate(draw = drawID),
paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw', drawID, '.csv'),
row.names = FALSE)
}
# Plot the first draw's errors:
plotting_tbl <- read.csv(paste0('generated_ensembles/', esm_name ,'/jump_errors_transitionyears_',esm_experiment,'_draw1.csv'),
stringsAsFactors = FALSE)
kable(plotting_tbl)
| tol | archive | draw | n_stitched | mean_gen | var_gen | mean_target | var_target | bias | E1 | E2 |
|---|---|---|---|---|---|---|---|---|---|---|
| 0.1 | with_ssp370 | 1 | 10 | 0.0392214 | 0.2218219 | 0.0147099 | 0.1754115 | 0.0245115 | 0.1397371 | 1.264580 |
| 0.1 | without_ssp370 | 1 | 30 | 0.0211943 | 0.2315472 | 0.0147099 | 0.1754115 | 0.0064844 | 0.0369666 | 1.320023 |
ggplot(data = plotting_tbl, aes(x = E1, y = E2, color = factor(tol), shape = archive)) +
geom_point() +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 0) +
ggtitle('Error metrics for single Draw 1, comparing window transition years of \ngenerated jumps to all years of target jumps')
Aggregated across draws:
# read in metrics for every draw
files <- list.files(path = paste0('generated_ensembles/', esm_name),
pattern = paste0('jump_errors_transitionyears_',esm_experiment ),
full.names = TRUE)
generated_metrics <- list()
for (i in 1:length(files)){
generated_metrics[[i]] <- read.csv(files[i], stringsAsFactors = FALSE)
}
generated_metrics <- do.call(rbind, generated_metrics)
generated_metrics %>%
group_by(tol, archive) %>%
summarize(aggregate_E1 = mean(E1),
aggregate_E2 = mean(E2)) %>%
ungroup %>%
mutate(total_draws = length(unique(generated_metrics$draw))) ->
aggregate_metrics
kable(aggregate_metrics)
| tol | archive | aggregate_E1 | aggregate_E2 | total_draws |
|---|---|---|---|---|
| 0.1 | with_ssp370 | 0.0891593 | 1.302006 | 10 |
| 0.1 | without_ssp370 | 0.0544899 | 1.292340 | 10 |
ggplot(data = aggregate_metrics, aes(x = aggregate_E1, y = aggregate_E2, color = factor(tol), shape = archive)) +
geom_point() +
geom_hline(yintercept = 1) +
geom_vline(xintercept = 0) +
ggtitle(paste('Aggregate errors comparing window transition years of \ngenerated jumps to all years of target jumps \naggregated across Ndraws ', unique(aggregate_metrics$total_draws) ))